{-
    Statistics.hs: Statistical functions.
    ----------------------------------------------------------------------------
    Version: 1.0alpha1
    Author:  Christopher Granade (cgranade@gmail.com)
    ----------------------------------------------------------------------------
    This file is a part of the Die Rolling Analysis Toolkit (DRAT) project.
    http://drat.googlecode.com/
    ----------------------------------------------------------------------------
    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 2 or 3 of the License.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
    ----------------------------------------------------------------------------
    (c) 2007 Christopher E. Granade
-}
module Statistics where
    import List

    describe    :: (RealFloat a) => [a] -> ((a, a, a), a, a, a)
    describe xs =  (quartiles xs, mean xs, std xs, skew xs)

    hist    :: (Ord a) => [a] -> [Int]
    hist xs =  (map length . group . sort) xs
    
    histNorm    :: (Ord a, Fractional b) => [a] -> [b]
    histNorm xs =  map (\x -> fromIntegral x / fromIntegral (sum (hist xs))) (hist xs)
    
    var    :: (Fractional a) => [a] -> a
    var xs = let n = fromIntegral $ length xs - 1 in
    	(1 / (n - 1)) * (sum (map (\x -> (x - mu)^^2) xs))
   	  where mu = mean xs
        -- (1 / (n - 1)) * (sum (map (^^2) xs)) - (n / (n - 1)) * (mean xs)^^2
    
    std    :: (Floating a) => [a] -> a 
    std xs = sqrt (var xs)
        
    rank :: (Fractional b, Ord b) => Int -> [b] -> b
    rank n xs = (sort xs)!!(n - 1)
        
    median :: (Fractional a, Ord a) => [a] -> a
    median xs
    	| (even . length) xs = (rank (halfn) xs + rank (halfn + 1) xs) / 2.0 
		| otherwise = rank (halfn + 1) xs
	  where
	    n = (length xs)
	    halfn = floor ((fromIntegral n) / 2)
    
    halflist :: [a] -> ([a], [a])
    halflist xs = (take halfn xs, drop (halfn + n `rem` 2) xs)
      where
        n = length xs
        halfn = floor ((fromIntegral n) / 2)
    
    quartiles :: (Fractional a, Ord a) => [a] -> (a, a, a)
    quartiles xs = (median lower, median xs, median upper)
      where
        (lower, upper) = halflist (sort xs)
    	
    
    mean    :: Fractional a => [a] -> a
    mean xs =  (sum xs) / (fromIntegral (length xs))
    
    {--
        I am not familiar enough with sample moments to guarantee that this is
        all correct.
    --}
    centralMoment      :: (Fractional a, Integral b) => [a] -> b -> a
    centralMoment xs k =  let mu = mean xs in
        (1 / (fromIntegral $ length xs)) * (sum (map (\x -> (x - mu)^^k) xs))
    
    stdCentralMoment      :: (RealFloat a, Integral b) => [a] -> b -> a    
    stdCentralMoment xs k =  (centralMoment xs k) / (std xs)^^(k)
    
    skew    :: (RealFloat a) => [a] -> a
    skew xs =  stdCentralMoment xs 3
    
    kurt    :: (RealFloat a) => [a] -> a
    kurt xs =  stdCentralMoment xs 4